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An overview of the capabilities of the CHarring Ablator Response (CHAR) code is presented. CHAR is a 
one-, two-, and three-dimensional unstructured continuous Galerkin finite-element heat conduction and ab- 
lation solver with both direct and inverse modes. Additonally, CHAR includes a coupled linear thermoelastic 
solver for determination of internal stresses induced from the temperature field and surface loading. Back- 
ground on the development process, governing equations, material models, discretization techniques, and nu- 
merical methods is provided. Special focus is put on the available boundary conditions including thermochem- 
ical ablation, surface-to-surface radiation exchange, and flowfield coupling. Finally, a discussion of ongoing 
development efforts is presented. 


Nomenclature 


ay, 34, temporal finite difference weights (sec~') 


B extent of reaction 

U vector of discrete unkowns 

v velocity (™/sec) 

w mass source (k&/m?-sec) 

m mass flux (kg/m?-sec) 

Q energy source (W/m?) 

q heat flux (W/m?) 

T volume fraction in virgin material 
n unit normal vector 

Lb dynamic viscosity (Pa - sec) 

@ porosity 

w basis function 

p density (k&/m?) 

K permeability tensor (m7) 

k thermal conductivity tensor (W/m-k) 
E activation energy (J/kg) or error 
€o total internal energy (J/kg) 

he total enthalpy (J/kg) 

k; pre-exponential factor (sec!) 
m reaction order 

ne number of components 

P pressure (Pa) 





* Applied Aeroscience and CFD Branch, Member AIAA. 
+ Applied Aeroscience and CFD Branch, Member AIAA. 


1 of 11 


American Institute of Aeronautics and Astronautics 


etyys 


temporal order of accuracy 
gas constant (J/kg-K) 
temperature (K) 

time (sec) 

test function 


Subscripts and Superscripts 


SS oat o 


quantity of the char 

quantity of the gas 

denotes finite dimensional approximation 
component index 

quantity of the mesh 

time index 

quantity of the virgin 


I. Introduction 


CHAR’s primary purpose is to perform engineering analysis for design work. Provide brief literature review. 
Provide snapshot of paper. 


II. Governing Equations 


CHAR solves the governing equations for heat conduction and gas flow through decomposing porous media in one, 
two, or three dimensions. In general, this would include 


e Continuity equations for each of the chemically reacting gaseous species resulting from flowfielnd injection and 


pyrolysis gas generation through in-depth decomposition 
Momentum equations accounting for porous media transport 
Energy equations for both the solid and gas phases 


Solid mass conservation equations for progression of pyrolysis 


However, CHAR employs several simplifying assumptions to reduce the equation set to be more amenable for engi- 
neering application and design analyses. These assumptions include 


Flowfield gases are not allowed to penetrate the surface. Consequently, all the gas inside the porouse medium is 
generated through pryolysis. 


Pyrolysis gases do not react with or condense on the porous material. 

Pyrolysis gas diffusion is negelected such that the elemental composition of the gas is constant. 
All ablation occurs at the surface, and volume ablation effects are neglected. 

The solid and gas are assumed to be at the same temperature. 

The gas is assumed to be a mixture of perfect gases in chemical equilibrium. 


The momentum equations can be sufficiently represented by the steady form of Darcy’s law. 


Under these assumptions, the governing equations for a moving mesh simplify to 








O(peo ~ : 
Energy: p60) —V.- (vr) + V - (dpgho,¥g) — Um V(pec) — Q = 0 
ald ’ node (1) 
Gas Mass: Ar — Wg + V-(bpg¥g) — Um: V(bpq) = 0 
node 
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where the node subscript denotes conservation with respect to a moving node and Darcy’s law defines the velocity 
field by 
Rvp 2) 
vg =-— 
7 oh 


Additionally, the solid mass conservation equation, given by 


Ops 
Ot 





ar (3) 


is solved on a fixed mesh and treated as a constitutive relation when defining the temporal derivative in the energy 
equation. 


Ill. Numerical Methods 


Il.A. Finite Element Formulation 


In order to develop the finite element equations, a Galerkin weak statement can be developed for the energy equation 
in Eqn. (1), by first multiplying it by a suitable test function, v, and integrating over the domain 22 while integrating 
the second and third terms by parts to give the natural boundary condition terms. The weak statement is then: Find 
pe, € H' such that 





ip fo teee + Vv- (vr) —Vv- (dpgho, Vg) —UUm > V (peo) — 10| dQ 


+ ¢ (uho, 1g + Vdcona,) I =0 Wwe Hy (4) 
r 
where the boundary mass flux due to gas convection is 


Mg = (pq) Ug: A (5) 


and the boundary heat flux is 
deond, =-kVT-n (6) 


Likewise, a Galerkin weak statement can be developed for the gas mass conservation equation in Eqn. (1). Again, 
the equation will be multiplied by a suitable test function, v, and integrated over the domain ( while integrating the 
third term by parts to give the natural boundary condition term. The weak statement is then: Find ¢p, € H' such that 





| (“ee — Vu- (Gpg¥g) — VUm > V(bpg) + inv) dQ + ¢ vmgdl =0 Vue HG (7) 
Q i 


Eqns. (4) and (7) can be discretized by expanding the independent variables and test functions in terms of a finite 


dimensional basis 
# nodes 


(p€0), (@) = S~ (peo), bj (a) (8) 


(dpa) -y (ps); (9) 


# nodes 


= 
a 
I 

M4 


Tj th; (x) (10) 


j=l 


# nodes 


Pi, ey= Pea) (11) 


j=1 


# nodes 


vn (#) = D> vidi (a) (12) 


i=l 
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where the subscript / is introduced to denote a finite dimensional approximation. The unknowns in the system are 
chosen to be the nodal temperatures and pressure, 7’; and P;. Since the unknowns are not functions of x, the partial 
differential equation system is reduced to an ordinary differential equation system in which the temporal derivatives 
can be defined as 








Mate AMC) 


dab, 
Ot é dt v5 (@) G2) 
j=l 
(Opa), "SS EGP a) 5 
OT a ag (14) 


Since the equation system should be satisfied for all combinations of nodal shape function coefficients, v;, their choice 
is arbitrary as long as a unique combination is chosen for each node. In typical finite element fashion, the coefficients 
are chosen to be v; = 5;, for the /*” nodal equation. Consequently, the nodal residual equations resulting from Eqns. (4) 
and (7) can now be written as 


[oe nans | Vos kV TAO — | ($Pq)p, ho, Vivi + VgdQ 
Q Q Q 
— [ vevm (peo), $ ideona, al + $ vihoyrngAl — f ¥.Qa0=0 5) 
Q T T Q 


and 


O(OPg)p, 


[gene _ i (pg), Vi + VgdQ — [bem -V(dbpq),d2 


+ Vangdt + | Wjw.dQ =0 (16) 
r 2 
for 7 = 1,2,...,# nodes. 


IlI.B. Time Discretization 


The semidiscrete weak form of the system given by Eqns. (15) and (16) is discretized in time using backward finite- 
difference schemes. Both first- and second-order accurate schemes may be derived from Taylor series expansions in 
time about U,,,,; as outlined by Kirk! to give 


O41 
ot 





= aU ys + BiUn + “WU n-1 + O (Ar (17) 


where the weights are given for p = 1 and p = 2 in Table 1. 


Table 1. First- and Second-Order Accurate Time Discretization Coefficients 











P ant Bt Vt 
1 -1 
Bia Bing p 
1 aco Atn41 
2 Be— Ve lat 11 | Atn Atn(Atn41+Atn) 











The resulting temporal derivatives are given by 





d (peo). i Fe fe 
— z = OF (peo); +4 Br (peo); + Yt (peo); ; (18) 
and 
d e } n n n- 
(909) da = a (dPq); a +B (dq); + Yt (Pq); 2 (19) 


which can be substituted into Eqns. (13) and (14) respectively. 
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IV. Material Model 


The material model is an extension of the CMA model developed by Moyer and Rindal.? It is assumed that all the 
pores are interconnected, and therefore pyrolysis gases occupy all of the pore space and are free to flow through it. 
Consequently, the density of the solid/gas system is described by 


p = $69 + Ps (20) 


where the solid density is a bulk density, the gas density is with respect to the space the gas occupies (pore space), and 
the porosity is equal to the gas volume fraction. In terms of units, Eqn. (20) can be expressed as 


¢ Pg Ps 


oO oOo nara stn OO 
[total mass] J [pore vol] [gas mass] _ [solid mass] 





= (21) 
[total vol] [total vol] [pore vol] [total vol] 


The solid material is modelled as a mixture of components each with their own rate law governing its decomposition. 
ne 
= Sol ips, (22) 


The units associated with the solid bulk density model in Eqn. (22) are 


Ps Ty Pi 





—— Canna ay ee th 
[solid mass] -¥ [initial vol of 2°" comp.] [mass of 7°” comp.] 


[total vol] C 


[total vol] [initial vol of qe comp.] 


It is assumed that the solid does not change volume due to thermal expansion, and therefore the total volume is 
constant. It is important to note that the solid description in Eqn. (22) is only a modeling assumption, and the solid 
is not truly comprised of nc components, species, or distinguishable materials. Taking the temporal derivative of 
Eqn. (22) gives the solid decomposition rate in terms of component decomposition rates. 


“hs =, = ye Pe ae -y Diy, (24) 


i= 





It is assumed that the decomposition of each component can be described by an Arrhenius relationship of the form 


psy. oe ea)" e—Hi/RT (25) 


which applies at a constant spatial location (as apposed to a given node that can move during the solution process). 
The decomposition kinetics in Eqn. (25) are simplified to ordinary differential equations and integrated as part of the 
solution of the energy equation. 

Since most thermophysical properties of the solid are only known for the virgin and fully charred states, the inter- 
mediate solid is modeled as some interpolated state between virgin and char. This interpolated state is characterized 
by the extent of reaction (3), or degree of char, given by 


B= Pu — Ps (26) 
Pu — Pe 


where the virgin and char bulk densities are known constants. The definition in Eqn. (26) can be rearranged to more 
clearly describe the interpolated state. 


Ps = (1 _ 6) Pu + Bic (27) 


Although the virgin and char materials are not distinguishable entities within the intermediate solid, Eqn. (27) reveals 
that the degree of char represents an effective char volume fraction within the solid (not in the solid/gas mixture). In 
a similar light, CMA defines an effective virgin mass fraction given by 


yy = (1 - ee) (28) 
Pu — Pe Ps 
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which can be related to the extent of reaction through 


yw = 2 (1-8) (29) 
Ps 
Similarly the char mass fraction is given by 
Ye=1-y == (30) 
Ps 


These effective parameters are used to interpolate between virgin and char properties for partially decomposed mate- 
rial. 


V. Code Architecture 


To expedite the software development process, the development team has chosen to take advantage of pre-existing, 
high quality, and freely available software libraries to avoid duplicating the efforts of others. Consequently, the CHAR 
program can be thought of as an ablation and thermal physics wrapper around a hierarchy of libraries that aid in the 
solution and data management processes. The primary software packages included in CHAR are outlined in Figure 1, 
and there are additional libraries that the noted libraries rely on for the functionality in CHAR. The majority of the 
physics and application-dependent code is in the top-level routines denoted by the CHAR box. This includes simulation 
management, initialization, calls to library functions, matrix assembly routines for element-level integration, boundary 
condition application, and material property data structures and look-up routines. Those routines comprised the bulk 
of the development effort for the current study. 

Eigen? is a library for linear algebra including: matrix/vector math, numerical solvers, and other related algorithms. 
It is primarily used in CHAR for its matrix and vector data structures and solutions to small linear systems 

Cantera? is a collection of object-oriented software tools for problems involving chemical kinetics, thermodynam- 
ics, and transport processes. For CHAR, it is used to calculate thermodynamic state properties and transport properties 
for gaseous mixtures in equilibrium. It is used for both the in-depth gas mixture and the environment gas mixture at 
surface conditions. Since repetitive calls to the surface equilibrium calculator can get expensive over the coarse of 
a simulation, Cantera is only used in the beginning of a simulation to build the required data tables over a range of 
relevant conditions for the problem. Subsequently, CHAR interpolates in these data tables to determine gas properties. 

The workhorse for CHAR is the 1ibMesh finite-element library, which is developed at NASA Johnson Space 
Center and the University of Texas at Austin with many other contributors around the world. 1ibMesh performs 
nearly all of the physics independent tasks of the software and can be thought of as the framework on which the 
physics-dependent application code is built. While the functions of 1ibMesh within CHAR are many, the notable 
functions include linear system data management, mesh data management, parallel communication, finite-element 
utilities that streamline integration, and adaptive mesh refinement. 1ibMesh also provides user-transparent interfaces 
to the Portable, Extensible Toolkit for Scientific Comuptation (PETSc) library,” which solves the assembled linear 
system, and METIS, ° which provides the domain decomposition for parallel execution. 


VI. Boundary Conditions 


VI.A. Thermochemical Ablation 
VI.B. Surface-to-Surface Radiation Exchange 


VI.C. Flowfield Coupling 
VII. Mesh Motion 


VIII. Adaptive Mesh Refinement 
IX. Inverse Solver 


X. Thermoelastic Solver 


CHAR has the ability to solve the linear thermoelastic equations in a one-way coupled quasi-steady fashion to give 
the displacement, stress, and strain fields for a body under thermoelastic loading. The equations are linearized by 
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Figure 1. CHAR software architecture. 
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only considering the linear terms in the strain-displacement relationships and by assuming that the material properties 
are not functions of the displacement. The thermoelastic equations are one-way coupled to the thermal and ablation 
equations such that an instantaneous temperature field provided from the thermal and ablation solver is assumed to be 
fixed during the thermoelastic solve. There is no feedback from the thermoelastic solver to the energy equation either 
in the form of solid deformation or thermoelastic coupling terms in the energy equation. The thermoelastic solver can 
be thought of as a post-processing of the thermal and ablation equations that can occur at a user-defined frequency. 
In the current implementation, the pressure field within a decomposing ablator is not considered to contribute to 
the thermoelastic response. It is recognized that this capability must be added to predict internal pressure induced 
spallation. 

The solution to the thermoelastic equations is quasi-steady in that the displacement field, uw, is assumed to be static 
at the end of each solve. This implies that 3; 

Uw 
aS 0 (31) 

The physical assumption is that the material elastically responds infinitely fast relative to changes in the temperature 
field. 


XI. Governing Equations 


The equations that govern the steady thermoelastic response of a solid under loading are 
Vo+f=0 (32) 


where a is the stress tensor, f is the body force per unit volume, p is the material’s density, and w is the displacement 
field. 
U1 
uU= luo (33) 


U3 


The symmetric stress tensor is defined as 


O11 O12 O13 


OF ; 022 023 =o 02 03 (34) 


033 


where the diagonal terms are the normal stresses and the off-diagonal terms are the shear stresses. where 


O1d 
Od = |o0a ford = 1,2,3 (35) 


03d 


The linear elasticity equation for each dimension can now be written as 
V-oat+fa=0 ford=1,2,3 (36) 


In order to derive the relationship between the stress and the displacements, it is necessary to recognize that the 
orthogonal 123-coordinate system that the problem is solved in (denoted as the “problem coordinate system” from here 
on) is in general not the same orthogonal 1’2’3’-coordinate system that the material properties are defined in (denoted 
as the “material coordinate system” from here on). Any value, vector, or tensor followed by the prime symbol, “’’, will 
be with respect to the material coordinate system, while any non-primed value, vector, or tensor will be with respect 
to the problem coordinate system. 

In the material coordinate system, the stress is related to the strain via Hooke’s law for small deformations 


o’ = K’': [e’ — ATa’] (37) 
where the strain, €’, is also a symmetric tensor and is defined via the linear strain-displacement relationship 
ie 


Bees ; [Val + (Vw (38) 
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ATa’ are the thermal strains where 


a, 0 0 
ai=|. a 0 (39) 
as 


and a’, are the coefficients of thermal expansion. K’ is the fourth-order stiffness tensor containing the 81 elastic 
coefficients, k rv Which are properties of the material and are defined with respect to the material’s coordinate system. 
Because of the linear assumption, the elastic coefficients are independent of the displacement, but they do vary with 
temperature. Due to symmetry, there are only 21 independent elastic coefficients in the stiffness tensor. Consequently, 
Voigt notation is often used to represent the stress and strain tensors to reduce the order of the problem. 


/ / / 
O11 11 ay 
/ / / 
022 €99 O92 
ot, é, al, 
o=| 8) =] 8] a’ = | 3 (40) 
093 2€53 0 
O44 265, 0 
O19 2€45 0 
and the stiffness tensor is simplified to 
/ / / / / / / / / / / / 
Ryva ky 122 ki133 ki 193 kyi31 ki Ki, Ky». K43 Big Kis Ki6 
/ / / / / / / 
kooi1 ko999 k5933 k5993 keo931 koo12 ‘ Ko. Ko3 Koa Ko5 Kae 
/ / / / / / / / / / 
K' = kegsiy k3399 3333 k3393 k3331 k3319 = . . 33 34 35 F36 (41) 
r / / / / / / / / / 
Ro311 k5399 k5333 k5393 k5331 ko319 ; ; ; B44 Fs Kg 
/ / / / / / / / 
3111 3122 3133 3123 3131 3112 ; 55 56 
/ / / / / z ; ‘ é ; é K! 
1211 1222 1233 1223 1231 1212 66 


where the subscript r is introduced to differentiate the reduced order tensor from the original 4°” order stiffness tensor, 
K’. With this new notation, the stress can be calculated according to 


o’ = K' |e’ — ATa’ | (42) 
or 
oy Ki, Kig Kis Kis Kis Ki6 ey aly 
Oho < Kho Ko3 Ko4 Kos Ke €59 al, 
033 ss : . K33 K34 K35 K36 €33 _—AT as (43) 
73 ; Ki, Kis Ke 2€43 0 
O13" : 55 56 2€31 0 
T2 : : : : : 66 2€4 0 


The stress vector, o’, can now be used to fill out the stress tensor, 0’. Next, a rotation matrix, C’, is required to get the 
stress tensor from the material coordinate system to the problem coordinate system. 


a =C'o'C (44) 


where C is defined in Section ??. 

In CHAR, the approach will be to use Eqns. (42) and (44) to define the stress tensor for integration of the PDEs. 
The strain vector in Eqn. (42) is defined using the displacements in the material coordinate system, u’. Since the 
unknowns are defined in the problem coordinate system, it is more convenient to work in the problem coordinate 
system when defining the strain vector. Consequently, the strain vector in Eqn. (42) will be calculated by rotating the 
strain tensor from the problem coordinate system to the material coordinate system. 


f <e =C™ ec"! (45) 


where the < symbol denotes filling in the strain vector given the values in the strain tensor. Since C' is an orthogonal 
matrix 


ct=sCc (46) 
and consequently 
é =e’ =CeC” (47) 
9 of 11 


American Institute of Aeronautics and Astronautics 


XII. Finite Element Formulation 


The Galerkin weak statement for the equations can be developed using Eqn. (36). The equation is multiplied by a 
suitable test function, v, and integrated over the domain, 22, while integrating the first term by parts to give the natural 
boundary condition term. The weak statement is then: Find u € H! such that 


Vo-oui— f ufatr— | v(oa-ayar =o Vu € Hy ford = 1,2,3 (48) 
Q Q r 


Eqn. (48) can be discretized by expanding the independent variables and test functions in terms of a finite dimensional 


basis 
# nodes 


(ua), (@) = S> (ua); by (x) (49) 


j=1 


# nodes 


V(ua), (@) = D> (wa); Vv; (w) (50) 


= 


# nodes 


vn (#) = S> vidi (a) (51) 


i=l 


# nodes 


Von (x) = So viVui (a) (52) 
i=1 
where the subscript / is introduced to denote a finite dimensional approximation. Now, the unknowns are no longer 
functions of xz. Since the equation system should be satisfied for all combinations of nodal test function coefficients, v;, 
their choice is arbitrary. Choosing a unique combination for each node allows for the development of a linear system 
where the number of equations equals the number of unknowns. In typical finite element fashion, the coefficients are 
chosen to be v; = 6,; for the /*” nodal equation. Consequently the nodal residual equations can now be written as 


Vii oadQ -f Wi fadQ — } Tydv =0 ford =1,2,3 (53) 
Q Q T 


for 7 = 1,2,...,# nodes. The surface traction is defined as 


Tij=0a-h ford=1,2,3 (54) 


XIII. Ongoing Development Efforts 


e Water phase change model 
e Unsteady porous flow momentum equations 
e Thermal non-equilibrium between solid and gas 


e Chemical non-equilibrium modeling 


XIV. Conclusions 
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